module evolve

  use dimensions

  complex, private      :: ExpKv(jz, jv)
  complex, private      :: ExpKx(jx, jDirection)
  !complex, private      :: ExpKVdt
  complex, private      :: ExpKz(jz)


contains

!-------------------------------------------------------------------------------------
subroutine evolve_init()
!-------------------------------------------------------------------------------------

  use parameters, only : dt, kz, vPhase, u, Dxxxx, &
                         iWaveFrame, iunit, isign

  implicit none

  integer              :: iDirection, ikz
  real                 :: v0


  !-- initialization --


  v0 = -vPhase

  ExpKz(:) = cexp( iunit * kZvec(:) * v0 * dt )  !  complex ExpKz(jz)
      
  do ikz=1,nkz
     ExpKv(ikz,:) = cexp( -iunit*vvec(:)*kZvec(ikz)*dt - Dxxxx * kZvec(ikz)**4 *dt    )
  enddo

  !  data iSign(1), iSign(2)/ 1, -1 /
  do iDirection = 1,2
     ExpKx(:,iDirection) = exp( - iunit * kXvec(:)*u*iSign(iDirection) * dt  )
  enddo

end subroutine evolve_init


!-------------------------------------------------------------------------------------
subroutine evolve_F_advect(dfu, &                              ! input/output
                           uxMat, f0vec)                       ! input
!-------------------------------------------------------------------------------------

  use fft_wrap, only   : fft_forward, fft_backward
  use parameters, only : Dxxxx, dt, iBcPerpPeriodic

  implicit none

  real,    intent(inout)   :: dfu(jz, jv, jx, jDirection)
  real,    intent(in)      :: uxMat(jx, jz, jDirection)
  real,    intent(in)      :: f0vec(jv)

  complex                  :: fz(jz)
  real                     :: aVar(jx), fVec(jx)
  integer                  :: ix, iz, iv, iDirection


 !------------------------------------------

    
  ! ExpKv(ikz,:) = cexp( -iunit*vvec(:)*kZvec(ikz)*dt - Dxxxx * kZvec(ikz)**4 *dt  )
   
  do iDirection = 1,2  

   
     !-- advect in z: df/dt = -v df/dz --

     do ix=1,jx
        do iv=1,jv
           fz(:) = dfu(:,iv,ix,iDirection)
           call fft_forward(jz,fz,fz)
           fz = fz/float(jz)
           call dealias(jz,fz) !  fz=FT(dfu)      
           fz(:) = fz(:)*ExpKv(:,iv)
           call fft_backward(jz,fz,fz)
           dfu(:,iv,ix,iDirection) = fz(:)
        enddo
     enddo
    

     !-- advect in x: d(f+f0)/dt = -d/dx [(f+f0) Ux(x)] --

     do iz=1,jz                      
        aVar(:) = uxMat(:,iz,iDirection)           
        do iv=1,jv !  complex fvc(jv),fvc0(jv)     fvc0(iv)=f0vec(iv)
           fVec(:) = dfu(iz,iv,:,iDirection) + 0.5 * f0vec(iv)   ! real fVec(jx)
           call UpWind2ndOrder_Variable_a( fVec,  dt, dx, aVar, jx, &
                0.5 * f0vec(iv),  0.5 * f0vec(iv), iBcPerpPeriodic )
           dfu(iz,iv,:,iDirection) = fVec(:)  - 0.5 * f0vec(iv) 
           !  0.5 f0vec  since each direction has 1/2 density,
           ! goes to 1/4 is 3D (1/6 with hexagonal symmetry)
        enddo
      enddo
     
   enddo ! iDirection

end subroutine evolve_F_advect

!-------------------------------------------------------------------------------------
subroutine evolve_F_accel(dfu, &                                ! input/output
                          Ez2D, fvc0T)                          ! input
!-------------------------------------------------------------------------------------

  use fft_wrap, only   : fft_forward, fft_backward
  use parameters, only : dt, Dvvvv, iunit

  implicit none

  real,    intent(inout)   :: dfu(jz, jv, jx, jDirection)
  real,    intent(in)      :: Ez2D(jx,jz)
  complex, intent(in)      :: fvc0T(jv)

  complex                  :: fvc(jv), z(jv)
  integer                  :: ix, iz, iv, iDirection


  !-- acceleration: df/dt = Ez df/dv

   do iDirection = 1,2 
      do ix=1,jx
         do iz=1,jz   !	begin spatial loop
            fvc(:) = dfu(iz,:,ix,iDirection)
            call fft_forward(jv,fvc,fvc)
            fvc = fvc/float(jv)
            call dealias(jv,fvc)
            z(:)   =  - Dvvvv*pvec(:)**4  - iunit*pvec(:)*Ez2D(ix,iz) + 1.E-09
            fvc(:) = fvc(:)*cexp(z(:)*dt) + &
                 0.5 * (1-cexp(z(:)*dt))*iunit*pvec(:)*Ez2D(ix,iz)*fvc0T(:)/z(:)
                 ! 0.5 since each direction has 1/2 density, 
                 ! goes to 1/4 is 3D (1/6 with hexagonal symmetry)
!	         do ip=1,np !  but not ip = 1 since pvec(1) = z(1) = 0 
!	         fvc(ip)=fvc(ip)*cexp(z(ip)*dt) + &
!                         (1-cexp(z(ip)*dt))*iunit*pvec(ip)*Ez2D(ix,iz)*fvc0T(ip)/z(ip)
!                end do
            call fft_backward(jv,fvc,fvc)
            dfu(iz,:,ix,iDirection)=fvc(:)       !   f(iv, iDirection) = dfu(iv, iDirection) + 0.5 * f0vec(iv)
         end do !		end  z spatial loop
      enddo       !       end  x spatial loop
   enddo        !       end idirection loop
   

end subroutine evolve_F_accel

!-------------------------------------------------------------------------------------
 
subroutine evolve_U(uxMat, uzMat, &                            ! input/output
                    dfu, Ex2D, &                               ! input
                    ExFilter, f0vec, uPerpMod)                 ! input

!-------------------------------------------------------------------------------------

  use fft_wrap, only   : fft_forward, fft_backward

  use parameters, only : u, dt, &
                         iWaveFrame, iBcPerpPeriodic, iFourier, &
                         pi, iunit, isign

  implicit none


  !-- input variables (from main code) --


  real,    intent(inout)   :: uxMat(jx, jz, jDirection)
  real,    intent(inout)   :: uzMat(jx, jz, jDirection)

  real,    intent(in)      :: dfu(jz, jv, jx, jDirection)
  real,    intent(in)      :: Ex2D(jx,jz)

  real,    intent(in)      :: ExFilter(jx)
  real,    intent(in)      :: f0vec(jv)
 
  complex, intent(in)      :: uPerpMod(jx)


  !-- local variables -----------------------

  complex    :: Phi1Dz(jz)
  complex    :: du(jx)

  real       :: fVec(jx)

  real       :: fVsVvec(jv)
  real       :: Phi1D(jx)
  real       :: uTemp(jz)
  real       :: uzVec(jz)
  real       :: uPerpVec(jz)

  integer    :: ix, iz, iv, iDirection



  !-- integrate over "v" to find longtitudinal velocity Uz -- 
       
  do ix = 1,jx
     do iz = 1,jz 
        do iDirection = 1,2
           !real dfu(jz,jv,jx,jDirection), fVec(jx), fVsVvec(jv),vvec(jv), uzMat(jx,jz,jDirection)
           fVsVvec(:) = dfu(iz,:,ix,iDirection) + 0.5 * f0vec(:)
           uzMat(ix,iz,iDirection) = sum( fVsVvec * vVec ) / sum( fVsVvec )
        enddo
     enddo
  enddo
   
   
  !-- advect u_perp in z:  dUx/dt = Uz dUx/dz + Ex -- 
    
  do iDirection = 1,2
   
     if(iWaveFrame==1)  then  !  advect ux in z since <u_z> approximately = -vphase =  V0, added on June 2, 2009
        do ix=1,jx  !  complex Phi1Dz(jz), ExpKz(:) = cexp( iunit * kZvec(:) * V0 * dt )
           if( iFourier==1 ) then    !  Fourier with linearization in uz

              Phi1Dz(:) = uxMat(ix,:,iDirection)  !  real uxMat(jx,jz,jDirection)
              call fft_forward(jz,Phi1Dz,Phi1Dz)
              Phi1Dz=Phi1Dz/float(jz)
              call dealias(jz,Phi1Dz) 
              Phi1Dz(:) = ExpKz(:) * Phi1Dz(:)
              call fft_backward(jz,Phi1Dz,Phi1Dz)
              uxMat(ix,:,iDirection) = Phi1Dz(:)
           endif
           if( iFourier.EQ.-1 ) then   !  upwind
              uzVec(:) =    uzMat(ix,:,iDirection)  !  advecting velocity in z direction, at each x
              uPerpVec(:) = uxMat(ix,:,iDirection)
              call SubUpwind_NonCon_Form( uPerpVec, dt, dz, uzVec, jz, 0.0,  0.0, 1, 1, uTemp ) !  periodic in z so far
              uxMat(ix,:,iDirection)= uPerpVec(:)
           endif
           if( iFourier.EQ.0 ) then   !  centered difference
              uzVec(:) =    uzMat(ix,:,iDirection)  !  advecting velocity in z direction, at each x
              uPerpVec(:) = uxMat(ix,:,iDirection)
              call SubUpwind_NonCon_Form( uPerpVec, dt, dz, uzVec, jz, 0.0,  0.0, 1, 0, uTemp ) !  periodic in z so far
              uxMat(ix,:,iDirection)= uPerpVec(:)
           endif
        enddo   ! real  uTemp(jz), uzVec(jz), uzMat(jx,jz,jDirection), , uPerpVec(jz)
     endif  ! else need to impliment centered difference variant of SubUpwind_NonCon_Form
     
     uxMat(:,:,iDirection) = uxMat(:,:,iDirection) + Ex2D(:,:) * dt  ! uxMat(ix,iz,iDirection)  all ix and all iz
         
  enddo        !       end idirection loop
 

  !-- transverse self-advection of u_perp: dUx/dt = -Ux dUx/dt  --       

  do iDirection = 1,2
      
     do iz=1,jz  !  real f0vec(jv)
        fVec(:) = uxMat(:,iz,iDirection)  !  real fVec(jx), uxMat(jx,jz,jDirection)
        call NonLinearAdvection( fVec, dt, dx, jx, u*iSign(iDirection),  u*iSign(iDirection), iBcPerpPeriodic )  
        !  MAYBE NOT periodic bc for ux
        
        du = fVec - u * iSign(iDirection) !  approx to uniform part of du, improve upon this,  complex du(jx)
        du = du * uPerpMod  !  make du periodic in x
        call fft_forward(jx,du,du)
        du=du/float(jx)
        call deAlias(jx,du)
        du = du * ExFilter
        call fft_backward(jx,du,du)
        fVec = du +  u * iSign(iDirection)                      
        uxMat(:,iz,iDirection) = fVec(:)
     enddo  ! iz
      
  enddo      ! iDirection loop
                              

end subroutine evolve_U
 



!-------------------------------------------------------------------------------------

end module evolve
